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1 Introduction Photonic crystals consist of a periodic arrangement of materials with different refrac- 
tive indices, usually a dielectric material and air with a period on the scale of the wavelength [I J. This 
structure allows to control light propagation on very short distances. Due to their small size photonic 
crystal structures can be used to fabricate integrated optical components like resonators and waveguides. 
Photonic crystal fibers ISO are an important application example which utilizes the light guiding princi- 
ples of photonic crystals. They are used e.g. in nonlinear applications H for high power transmission of 
light in optical astronomy. 

Radiation losses from PhC fibers are one of the loss mechanism which potentially restrict their tech- 
nological application |5 |. They have therefore to be controlled by proper component design. Due to the 
complex geometry and large number of glass air interfaces the accurate computation of light propagation 
in photonic crystal structures is a challenging task Q. To achieve a realistic model of a real device its finite 
size and the exterior domain have to be taken into account. The coupling of a photonic crystal fiber to the 
exterior will always lead to unwanted losses which with a proper fiber design are very small. For numer- 
ical computation transparent boundary conditions are used to take into account the surrounding material. 
Even for simple geometries the accurate computation of radiation losses can become difficult and many 
numerical methods can fail especially for systems where these losses are very small \JJ. 

In this paper we review an approach to this problem using adaptive finite element algorithms and we 
apply this approach to accurate computation of radiation losses in hexgonal and kagome-structured PhC 
fibers. Computational results agree well with experimental transmission spectra measured in |5|. In the 
context of topics related to the DFG priority programme photonic crystals, we have applied adaptive finite 
element algorithms to a variety of nano-optical problems ll8ll9l [l^[T11[l^[T3l[T4l l6l[T51. 



2 Formulation of propagation mode problem For the analysis of photonic crystal fibers we start 
with the derivation of the mathematical formulation of the propagation mode problem of the electric and 
magnetic field. The geometry of a waveguide system like a photonic crystal fiber is invariant in one spatial 
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dimension along the fiber. Here we choose the z-direction. Then a propagating mode is a solution to 
the time harmonic Maxwell's equations with frequency uj, which exhibits a harmonic dependency in z- 
direction: 



E = Epin(x,?/)exp(?fc^2;) 
H ^ Hpin(a;,?/)exp(?fc^z) 



(1) 



Epm(2:, y) and Hpin(a;, y) are the electric and magnetic propagation modes and the parameter is called 
propagation constant. If the permittivity e and permeability /i can be written as: 
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we can split the propagation mode into a transversal and a longitudinal component: 
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Inserting ([T]) with Q and (O into Maxwell's equations yields: 
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Now we define Ez = kzEz and get: 



with 



E± 

E,. 



A = 
B = 



kiB 



El 

E, 



PV^Ai-iVi -P-c^^g^i 




Pt,l\P 







iW± ■ P/iJ^^P 



-iPnl\PVi_ 
• Pfil\PV± - ujhz 



(6) 

(7) 
(8) 



Eq. (l6|l is a generaUzed eigenvalue problem for the propagation constant kz and propagation mode Epm (a; , y ) . 
We get a similar equation for the magnetic field Hpin(a:, y) exchanging e and /i. For our numerical analysis 
we define the effective refractive index nofr which we will also refer to as eigenvalue: 
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where Aq is the vacuum wavelength of light. 
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3 Discretization of Maxwell's equations with the Finite Element Method For the numerical solution 
of the propagation mode problem Eq. (|6]l derived in the previous section we use the finite element method 
ri6J which we sketch briefly in the following. We start with the curl curl equation for the electric field. 
Since we want to solve an eigenvalue problem we are looking for pairs E and kz such that: 



1, 



-E = inn 



Vfc^ X -Vfc. X E - 

— Vfc xE)xn = F given on F (Neumann boundary condition) 
M ' J 



(10) 



(11) 



holds, with Vfc^ = [dx, dy,ikz]'^ . For application of the finite element method we have to derive a weak 
formulation of this equation. Therefore we multiply ( fTOb with a vectorial test function ^ ^ V = H{curl) 
lfT6l and integrate over the domain Vl: 



Vfe, X -Vfe^ X E 



* • E = , V* G 



(12) 



where bar denotes complex conjugation. After a partial integration we arrive at the weak formulation of 

Maxwell's equations: 

Find E e y = H{curl) such that 



(Vfc, X $) 



1, 



Vfc, X E 



in V \M 
We define the following bilinear functionals: 
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Now we discretize this equation by restricting the space V to a finite dimensional subspace Vh C V, 
(a) (b) (c) 
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Fig. 1 (a) Computational domain and (b) triangulation of hollow-core photonic crystal fiber; (c) example of vectorial 
ansatz function on a reference patch. 



dimVh = N . This subspace and therewith the approximate solution are constructed as follows. One starts 
with a computational domain Q,, see Fig. [11 a)- This domain is subdivided into small patches, e.g. triangles 
in 2D and tetrahedrons in 3D, FiglHb). On these patches vectorial ansatz functions Vi are defined. Usually 
on each patch the ansatz funtions Vi form a basis of a polynomial function space of a certain degree p 
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lfT6l . The approximate solution for the electric field is a superposition of these ansatz functions of all 
patches: 



N 



E/i = ^ aii^, (16) 
Together with (fT4l) . dTSI ) the discrete version of Maxwell's equations ( [13]) reads: 

N 

ci^a{v,,v,) = f{vj) , Vj = 1, . . . , iV (17) 



which is a linear system of equations for the unknown coefficients a^: 

( \ 

with ^ a{vi,Vj), fj = J{vj), a= \ ... (18) 

V ajv / 

The matrix entries a{vi, Vj) arise from computing integrals (fT4b . In practice these integrals are evaluated 
on a reference (unit) patch. Such a reference patch together with a vectorial ansatz function is shown in 

Fig.mc). 

In the above sketch we assumed for simplicity that boundary conditions ( fTTT ) are known for the electric 
field E. However here we want to take the infinite exterior into account. Therefore the eigenvalue problem 
^ has to be solved on an unbounded domain M?. This leads to the computation of leaky modes which 
enable us to estimate radiation losses. Since our computational domain still has to be of finite size, we apply 
so-called transparent boundary conditions to dO,. We realize these boundary conditions with the perfectly 
matched layer (PML) method [ 17 1. Details about our numerical implementation are described in 1 18|. The 
propagation constant (and the effective refractive index) becomes complex and the corresponding mode 
is damped according to exp(— while propagating along the fiber, see Eq. ([T]). 

Applying the finite element method to propagating mode computation has several advantages lfT9ll20l . 
The flexibility of triangulations allows the computation of virtually arbitrary structures without simplifi- 
cations or approximations, as illustrated in Fig. [Tfb) and|7jb). By choosing appropriate ansatz functions 
Vi{x, y) for the solution of Maxwell's equations, physical properties of the electric field like discontinu- 
ities or singularities can be modeled very accurately and do not give rise to numerical problems. Such 
discontinuities often appear at glass/air interfaces of photonic crystal fibers, see Fig. |2] Adaptive mesh- 
refinement strategies lead to very accurate results and small computational times. Furthermore the FEM 
approach converges with a fixed convergence rate towards the exact solution of Maxwell-type problems 
for decreasing mesh width (i.e. increasing number of sub-patches) of the triangulation. Therefore, it is 
easy to check if numerical results can be trusted 1 16 |. 

Especially for complicated geometrical structures the finite element method is better suited for mode 
computation than other methods. In contrast to the plane wave expansion (PWE) method, whose ansatz 
functions are spread over the whole computational domain (plane waves) the FEM method uses localized 
ansatz functions, see Fig. [Uc). In order to expand a solution with discontinuities as shown in Fig. |2]a large 
number of plane waves would be necessary using the PWE method. This leads to slow convergence and 
large computational times 1 12). 

For accurate and fast computation of leaky eigenmodes we have implemented several features into the 
FEM package JCMsuite. As we will see later e.g. high-order edge elements, a-posteriori error control and 
adaptive and goal-oriented mesh refinement increase the convergence of numerical results dramatically. 



4 Computation of leaky modes in hollow core photonic crystal fibers The only simpUfication we 
make for the computation of leaky modes is to extend the fiber cladding to infinity and thereby neglect its 



5 





finite size. This is justified if the cladding of the fiber is much larger than the microstructured core, and 
no light entering the cladding is reflected back into the core, which is usually the case. We investigate 
two different types of hollow-core photonic crystal fibers (HCPCF). The first type has a hollow core which 
corresponds to 19 omitted hexagonal cladding cells Fig. |3la) the second type is kagome-structured [SJ, Fig. 
Ob). The cross sections of the HCPCFs depicted in Fig. [3]have a Cgv invariance. Therefore it is possible 

(a) (b) 




Fig. 3 (a) Geometry of 19-cell HCPCF and (b) kagome-structured fiber used for mode computation 



to take only a fraction of the layout as computational domain. At the artificial inner boundaries appropriate 
boundary conditions have to be stated, setting the tangential magnetic and electric field to respectively. 
Table [T] shows the first and second eigenvalue for the full half and quarter fiber used as computational 
domain. The geometrical parameters are given in Fig. |2] The computed values are identical up to the 
chosen accuracy. The number of unknowns and computational time is of course much smaller for half and 
quarter fiber cross section. 

The imaginary and real parts of the eigenvalues given in Table[T]differ by up to 11 orders of magnitude. 
The complex eigenvalue leads to a dampening of the mode according to 

|E|2 cx e-^^C^^)^, (19) 

see Eq. ([T]l. Computing the fiber without transparent boundary conditions but setting the tangential com- 
ponent of the electric field to at the outer boundary we get nig = 0.99826580015 for the first eigenvalue, 
i.e. a real eigenvalue. Comparing to Table [T] we see that taking into account the finite size of the fiber 
and coupling to the exterior does not change the real part. However when analyzing radiation losses of a 
fiber, the imaginary part of the effective refractive index rics is the quantity of interest. Computation of 
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unknowns 1 st eigenvalue 2nd eigenvalue 

full fiber 3070641 0.99826580337+8.9297 • lO'^'H 0.992141777+ 1.3312 • 10-^"^ 

halffiber 1567377 0.99826580252 + 8.9311 • lO-^^i 0.992141758 + 1.3314 • 10^1°^ 

quarterfiber 781936 0.99826580254+ 8.9311 • lO^^^i 0.992141754+ 1.3315 • 10^1°^ 

Table 1 First, second and third eigenvalue computed with full, half and quarter fiber as computational domain. 



leaky modes with only small losses is therefore a multi-scale problem which is numerically very difficult 
to handle. In the next section we will introduce a goal-oriented error estimator which controls the grid 
refinement in order to obtain an accurate imaginary part as fast as possible. 



5 Goal oriented error estimator The finite element method enables us to refine patches of the un- 
structured grid only locally. We use error estimators to control the refinement process of the grid. Usually 
such an error estimator is defined by minimizing a target functional j(E), i.e. it is goal-oriented. The 
target functional depends on the solution of the electric field E. Since we are interested in radiation losses 
this will be the imaginary part of the propagation constant ifTSll . From Maxwell's equations applied to 
a waveguide structure one can derive the following expression for the imaginary part of the propagation 
constant: 



where 



Pn(E) 



2Po(E)' 
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2K(w) 
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E X -Vfc E 
E X -VkE 



(20) 



(21) 



(22) 



is the power flux of the electric field through the cross section of the computational domain and the in 
plane power flux across the boundary dil of the computational domain respectively. Equation (|20l i also 
shows us that the imaginary part of reflects radiation leakage from the photonic crystal fiber to the 
exterior. The right hand side of ( l20b is the target nonlinear functional j(E) which is used for the control of 
the error estimator The finite element mesh should now be adapted such that j(Eh) — i(E) is minimized, 
where E/i is the finite element solution and E the exact solution. This task can be embedded into the 
framework of optimal control theory lIZTl . We define the trivial optimization problem: 



j(E)-j(E,,) 



min 

*e H(curl) 



•j(E„) : a(*,*) = /(*)V*eH(cMr/)}, 



(23) 



This formulation is trivial because the restriction simply states that 5* is a solution of Maxwell's equation. 
The minima of ( |23]) correspond to stationary points of the Lagrangian density 

/:(E, E*) = j(E) - j(E,0 + /(E*) - a (E*, E) , 

where E* denotes the 'dual' variable (Lagrangian multiplier). Hence, we seek the solution (E, E*) to the 
Euler-Lagrange system 



a(*,E) = /(*) y^eH{curl) 
a(E*,*) = /(E;*) e H{curl). 



(24) 
(25) 



Equation (l24b is the variational form of the original Maxwell's equations. In the dual Eq. dZSl) the target 
functional j(E) = appears in form of its linearization j'(E). A finite element discretization of the 
Euler-Lagrange system yields a supplemental discrete problem 
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To quantify the eiTor of the finite element solution Eh we introduce the primal and dual residuals: 

p(Ehr) = f{-)-ai;Eh) (26) 
p{Elr) = fiEh:-)-a{El-). (27) 

The residuals quantify the error inserting the approximate solution into the exact non-discretized Maxwell's 
equations. For the exact solution one finds p(E;i; E) — Oandp(EJj;E) = 0. These residuals are computed 
for each patch and only patches with the largest residuals are refined. More details about the mathematical 
formulation and implementation can be found in ifTSl . A closer look to p{Eh; •) which is also used for 
field-energy based adaptive refinement will be given in the next section. 



6 Convergence of eigenvalues using different error estimators In this section we will look at the 
convergence of the computed eigenvalues. As computational domain we use the HCPCF depicted in Fig. 
[3 a). We compare three different refinement strategies namely uniform refinement and two different adap- 
tive goal-oriented refinement strategies. In an uniform refinement step each triangle is subdivided into 4 
smaller ones. In an adaptive refinement step a predefined resdiuum p of the numerical solution is mini- 
mized by refining only triangles with largest values for this residuum. The two different adaptive strategies 
we are using differ by the chosen residua. The first one only uses ( |26] |. We will refer to it as strategy A. 
The second strategy B which was introduced in the previous section uses (l27b in addition to ( |26] |. 

Let us have a closer look at the common residuum ( l26b of both strategies. We start with the weak 
formulation of Maxwell's equations (fTSl) on a sub patch Hi of our triangulation and reverse the partial 
integration: 







X — Vfc X E 



E 



-Vfe, X E 



X n 



d^r, (28) 



where Ti is the boundary of the sub patch fJ, and 



X E 



is the difference of the electric field 
E and the permeability p, on both sides of this boundary. Since Eq. ( [28] l holds for arbitrary $ the terms 
in brackets in both integrals have to vanish. For an exact solution to Maxwell's equations the first term 
vanishes because this is the Maxwell equation itself and the second term because the tangential component 
of ^ Vfe^ X E is continuous across a boundary. Since we approximate the exact solution by Eh both terms 
will generally not vanish. Therefore we define the residuum pi of sub patch ilf. 



Pi{Eh;Eh) ^ h^. 



Vfe, X -Vfc^ 



X Eh 



-Ey 



d^r+h. 



-V. X Eh 



X n 



d^r, 
(29) 



where hi is the size of the sub patch. For the solution Eh = E of Maxwell's equations we find pi (Eh ; E) = 
0. The residuum is therefore a measure how well the approximation Eh fulfills the exact Maxwell's equa- 
tions. 

Figure |4] shows the relative error of the fundamental eigenvalue nlff in dependence on the number of 
unknowns of the FEM computation for different refinement strategies and finite element degrees p. The real 
part of the effective refractive index converges for all finite element degrees and all refinement strategies, 
see Fig. |4](a),(c),(e). For higher finite element degrees we find faster convergence. For the lowest finite 
element order and coarsest grid (with « 16000 unknowns) we have a relative error of « 10^^. This error 
decreases down to 10^^^ for fourth order elements with 4- 10® unknowns. Figure|4](b),(d),(f) show that the 
imaginary part converges much slower. For uniform refinement and finite element degrees of p — 1, 2 we 
do not find convergence at all. The adaptive strategies A and B refinement also lead to poor convergence 
for low finite element degrees. For an accurate computation of small losses therefore high finite element 
degrees are necessary. For a relative error of lO^'^ = 0.1% we already need « 10® unknowns for p = 3, 4. 
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Fig. 4 Relative error of fundamental eigenvalue in dependence on number of unknowns of FEM computation for 
different refinement strategies and finite element degrees p. Parameters: A = 1550 nm, r = 300 nm, w = 50 nm, 
t = 170 nm, 6 cladding rings, wavelength A — 589 nm. Adaptive refinement strategy A minimizes residuum J26t , 
strategy B minimizes residua J26b and J27b 
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Fig. 5 Relative error of fundamental eigenvalue in dependence on number of unknowns of FEM computation for 
different refinement strategies and finite element degree p = 3. Parameters: A — 1550 nm, r — 300 nm, w — 50 nm, 
t = 170 nm, 6 cladding rings, wavelength A — 589 nm. (a) uniform refinement, (b) adaptive refinement, (c) goal- 
oriented refinement 



In Fig. |5ja) and (b) the different refinement strategies are compared for fixed finite element order p = A. 
For the real part we find fastest convergence for adaptive strategy A. As explained before the corresponding 
error estimator refines those triangles where the electric field has a large deviation from the exact solution. 
The error estimator of strategy B also uses this residuum but furthermore the residuum of the dual problem, 
Eq. ( |26l ) and ( |27] i. Therefore it also converges faster than uniform refinement but slower than strategy A. 
For the imaginary part the adaptive refinement strategy B shows fastest convergence. Strategy A shows 
almost no benefit for the first two refinement steps in comparison to uniform refinement. For the exact 
result for the convergence plots we used the most accurate result (p = 4) for the real and imaginary part 
obtained from strategy A and B respectively. 

Fig. |6] shows the benefit of high order finite elements. Here the adaptive refinement strategy A was 
used corresponding to Fig. |4tc),(d) but with finite element degree up to order p — 7. While the real part 
of the eigenvalue converges very fast even with order p = 2 the imaginary part can be computed much 
more accurate with higher order finite elements. Compared to p — 2 using orders greater than p = 4 the 
relative error is 2 orders of magnitude smaller for the same number of unknowns. Computational times for 
an 2.6 GHz AMD Opteron processor system are also shown. 

A comparison of the convergence and computation efficiency between our finite element package and 
the MIT Photonic-Bands (MPB) package (plane wave expansion method) is presented in ||22J where bloch 
modes of photonic crystal structures were computed. The FEM computations showed a much higher 
convergence rate. Results of the same accuracy could be computed over 100 times faster then with the 
MPB package. In [12] FEM computations of guided modes in HCPCFs are also compared to plane wave 
expansion (PWE) computations. Eigenmodes which were computed in about a minute with our finite 
element package took several hours with the PWE method. 

7 Optimization of HCPCF design Since we are enabled to compute radiation losses very accurately 
we can use our method to optimize the design of photonic crystal fibers to reduce radiation losses. The 
basic fiber layout is a 19-cell core with rings of hexagonal cladding cells. Fig. ^a). Since these cladding 
rings prevent leakage of radiation to the exterior for core guided modes we expect that an increasing 
number of cladding rings reduces radiation leakage and therefore reduces 3(ncff). This is confirmed by 
our numerical simulations shown in Fig. [8] The radiation leakage decreases exponentially with the number 
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(a) (b) 




Fig. 6 Relative error of fundamental eigenvalue in dependence on number of unknowns of FEM for adaptive refine- 
ment and different finite element degrees p. Parameters: A — 1550 nm, r = 300 nm, w = 50 nm, t = 170 nm, 6 
cladding rings, wavelength A — 600 nm. Computational times are shown. 



of cladding rings and thereby the thickness of the photonic crystal structure. This behavior agrees with 
the exponential dampening of light propagating through a photonic crystal structure with frequency in the 
photonic band gap. For our further analysis we fix the number of cladding rings to 6. The free geometrical 

(a) (b) 




Fig. 7 (a) geometrical parameters describing HCPCF: pitch A, hole edge radius r, strut thickness w, core surround 
thickness t; (b) detail from a triangulation of HCPCF. Due to the flexibility of triangulations all geometrical features 
of the HCPCF are resolved. 



parameters are the pitch A, hole edge radius r, strut thickness w, and core surround thickness t depicted in 
Fig- Eta) together with the triangulation |7tb). Fig. [8] shows the imaginary part of the effective refractive 
index in dependence on these parameters. For each scan all but one parameter were fixed. For the strut 
thickness w and the hole edge radius r we find well-defined optimal values which minimize 3(ncff)- For 
pitch A and core suiTound thickness t a large number of local minima and maxima can be seen. Now 
we want to optimize the fiber design using multidimensional optimization with the Nelder-Mead simplex 
method lfT2l . To reduce the number of optimization parameters we fix the hole edge radius to the determined 
minimum at r = 354 nm since its variation has the weakest effect on 3(ricff). For optimization we have 
to choose starting values for A, t and w. Since the simplex method searches for local minima we have to 
decide in which local minimum of A and t we want to search. We choose t = 152 nm since here 3(neff) 
has a global minimum and A — 1550 nm since the bandwidth of this minimum is much larger than for 
the global minimum at A = 1700nm. Optimization yields a minimum value of 3(riofi) = 5 • 10^^^^ 
for the imaginary part of the effective refractive index. The corresponding geometrical parameters are 
A = 1597 nm, w = 38 nm ,t = 151 nm. 
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(a) (b) 




r 

Fig. 8 Imaginary part of effective refractive index 3(ncft ) in dependence on: (a) number of cladding rings, (b) pitch 
A, (c) core surround thickness t, (d) strut thickness w, (e) hole edge radius r. Parameters: A = 1550 nm, r = 300 nm, 
w — 50 nm, t = 170 nm, 6 cladding rings, wavelength A — 589 nm. 

8 Kagome-structured fibers In iH attenuation spectra of large-pitch kagome-structured fibers have 
been measured experimentally. Fig. |9ja), (b) show the first and fourth fundamental core mode of such a 
19-cell kagome-structured fiber. The corresponding layout is given in Fig|9{c). In this section we compute 
attenuation spectra numerically. Since we only take into account radiation losses we do not expect quanti- 
tative agreement with experimental measurements. Discussions about other loss mechanism can be found 
in [12, 5|. We fix the geometrical parameters of the fiber and search for leaky eigenmodes in a wavelength 
interval. In the experiment core modes could be excited selectively. Here we also use the imaginary part 
of the fundamental core mode for the attenuation spectra. A small imaginary part then corresponds to low 
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(a) (b) (c) 




Fig. 9 (a) Fundamental and (b) fourtii core mode of 19-cell kagome-structured HCPCF (c). Parameters according to 
Tabled 

Layout Hollow-core Pitch A [/iin] Strut Width [fim] 
~K TO^celi 109 051 

B 1-cell 11.8 067 

Table 2 Layouts of kagome-structured fibers used for mode computation. 



losses and therefore high transmission. Furthermore we look at the confinement of the computed modes. 
Therefore we compute the energy flux of the mode within Score and outside i?strut the hollow core. A well 
confined mode then has a confinement close to 1. We expect that well confined modes have small 

in n I— i- 

losses. For simulation we use layouts corresponding to pi shown in Fig. Qlh) and|9lc) with parameters 
given in Table |2] 

(a) (b) 




0,6 0,8 1 1,2 1,4 "%,71 0,72 0,73 0,74 

A [^m] A [^m] 

Fig. 10 (a) Real part of effective refractive index of fundamental leaky mode in dependence on wavelengtli A for 
19-cell and 1-cell kagome-structured fiber. Parameters: see Table|2]- (b) zoom into Fig. I Uf a). 

Fig. [TOl' a) shows the real part of the effective refractive index of the fundamental mode in dependence 
on the wavelength. Since the core is filled with air it is close to 1 and hardly changes. Fig. [TT] shows the 
imaginary part of the eigenvalues which is proportional to the losses according to Eq. ( fT9] l. For the 19-cell 
fiber we find a region of low transmission in the wavelength interval 950nm-1050nm. In |5| a region of 
low transmission spans from 850nm-1050nm. For the 1-cell fiber we find a peak of high attenuation at 
A = 700 nm and high attenuation from 1200nm-1400nm. A dip in the transmission spectrum can also be 
found in the corresponding experimental measurement at 630nm. The range of the low transmission band 
however differs from the simulated. Experimentally it reaches from 800nm-1250nm. An important loss 
mechanism is coupling of the fundamental mode to interface modes at the glass air interfaces of the fiber 
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|l5][T2l. This is not taken into account in our simulations and could explain the mentioned disagreement. 
The regions of very high attenuation correspond to poorly confinement modes shown in Fig. [13] 

We notice that both attenuation spectra are very noisy. Since our computed results have converged we 
assume that this is no numerical artifact. To further investigate the large number of local extrema we zoom 
into the 19-cell spectrum from 0.7 15nm-0.735nm where a very large local peak can be seen, see Fig. [Tol' b). 
To explain the resonance in the attenuation spectrum at 725nm we look at the field distribution within the 
kagome structure. 

(a) (b) 




A [/im] A [^m] 

Fig. 11 Imaginary part of effective refractive index of fundamental leaky mode in dependence on wavelength A for 
(a) 19-cell and (b) 1-cell kagome- structured fiber. Parameters: see Table[2] 



Fig. [T2I shows the fundamental eigenmode for A — 719.25 nm where we find very low attenuation 
and for A = 725 nm at the local peak attenuation. The intensity of the eigenmodes in Fig. \T7\ a) and (b) 
is shown for the same pseudo color range. The field distribution for A = 725 nm is more intense in the 
glass struts which connect the core and the cladding of the fiber. Light from the fundamental core mode 
is coupled much stronger into the first neighbouring triangles of the kagome structure for A — 725 nm. 
They could be seen as resonators coupled to the hollow core and being excited by the core mode. Better 




Fig. 12 Pseudo color image of intensity of fundamental core mode for (a) A — 725 nm and (b) A — 719.25 nm. 
Parameters: Layout A, Table|2l 



14 



J. Pomplun, L. Zschiedrich, R. Klose, F. Schmidt, and S. Burger: FEM simulation of radiation losses 
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Fig. 13 (a) Fraction of electric field intensity of fundamental leaky mode located inside hollow core in dependence 
on wavelength A for (a) 19-cell and (b) 1-cell kagome-structured fiber. Parameters: see Table|2] 



coupling then leads to more light leaving the core and therefore higher attenuation. The coupling into the 
complicated kagome structure depends very sensitively on the wavelength and the shape of the mode. 

Finally Fig. \T3\a) shows the confinement of the fundamental mode in dependence on the wavelength. 
Regions of low loss correspond to regions of high ^{ucs), compare Fig. [TT] 

9 Conclusion We have investigated radiation losses in photonic crystal fibers. We have shown that 
with our FEM analysis we can efficiently determine the complex eigenvalues to a very high accuracy. The 
finite size of the photonic crystal fibers was taken into account by transparent boundary conditions reaUzed 
with the PML method. 

Computing leaky propagation modes in HCPCFs eigenvalues with real and imaginary part differing by 
over 10 orders of magnitude were found. The relative error of the imaginary part was thereby about 7 orders 
of magnitude larger than the relative eiTor of the real part. In order to precisely compute the much smaller 
imaginary part of the eigenvalues special techniques were implemented and applied to 2 different types of 
HCPCFs, namely HCPCFs with hexagonal cladding cells and kagome-structured fibers. A goal-oriented 
error estimator was introduced which focused on the accurate computation of a target functional, in our 
case the imaginary part of the eigenvalue. After each computation on a refinement level this error estimator 
was used to adaptively refine the grid. In a convergence analysis it was shown that the goal-oriented error 
estimator led to faster convergence of the quantity of interest compared to uniform refinement or standard 
field energy based refinement strategies. Also the usage of high order finite elements significantly increased 
the accuracy of the imaginary part. 

Since the imaginary part of a fiber could be computed in about 10 minutes with a relative error smaller 
than 10^'' it was possible to use the finite element method to automatically optimize a fiber design with 
respect to pitch, strut thickness, cladding meniscus radius and core surround thickness in order to minimize 
radiation losses. Furthermore attenuation spectra of a 1-cell and a 19-cell kagome-structured fiber were 
computed and compared to experimental results. It was shown that high losses were connected to poorly 
confined modes within the hollow-core. Furthermore the appearance of a large number of local minima 
and maxima in the attenuation spectra was explained by analyzing the intensity distribution of the core 
modes. 
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